Kolmogorov Turbulence in a Random Force 
Driven Burgers Equation 



Alexei Chekhlov and Victor Yakhot 

Program in Applied and Computational Mathematics, 
Princeton University, 
Princeton, New Jersey 08544, USA 
February 5, 2008 

Abstract 

The dynamics of velocity fluctuations, governed by the one-dimensional Burg- 
ers equation, driven by a white-in-time random force / with the spacial spectrum 



|/(fc)| oc is considered. High-resolution numerical experiments conducted in this 
work give the energy spectrum E{k) oc with (3 = | ± 0.02. The observed two-point 
correlation function C(k,uS) reveals u> oc k z with the "dynamic exponent" z ~ 2/3. 
High-order moments of velocity differences show strong intermittency and are domi- 
nated by powerfull large-scale shocks. The results are compared with predictions of 
the one-loop renormalized perturbation expansion. 

PACS number(s): 47.27.Gs. 



From a theoretical viewpoint, one of the most challenging features of strong hydrodynamic 
turbulence is an interplay between random almost gaussian background and coherent ordered 
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structures responsible for deviations from gaussian statistics. Although coherent structures 
were visualized in three-dimensional flows as sheets or tubes of high vorticity [|T|, little is 
known about their analytic structure, stability and, as a consequence, about their relevance 
to turbulence dynamics. In two-dimensional systems, the role of coherent structures is much 
better understood: the flow can be decomposed into two components, the background field 
having close to gaussian statistics and coherent, extremely stable point vortices, responsible 
for strongly non-gaussian features of the flow . Still, the analytic structure of the vortices 
and the distributions of their sizes and strengths are not yet understood and this is why the 
full statistical theory of two-dimensional turbulence does not exist. 
The analytic properties of the one-dimensional Burgers equation 

dv ldv 2 d 2 v 

subject to initial and boundary conditions, are understood rather well: the flow is dom- 
inated by the shocks, leading to E(k) oc k~ 2 ||. Moreover, in some cases, the Burgers 
equation has a stationary solution. For example, if v — —U and v — U at x — oo and 
— oo respectively, U(x) = —U tanh {x U/(2 vq)), which describes a single shock of the width 
I ~ u/U. In this particular solution, "fluid" particles, created at the boundaries, are carried 
towards the center of the shock where they disappear. The shock formation is the most 
significant dynamic property of the Burgers equation and the shocks were studied in the 
systems decaying from some initial conditions and driven by the large-scale random noise 
H. In the latter case, the energy spectrum is E(k) oc k~ 2 and all velocity structure functions 



S2n{r) = (u(x + r) — u(x)) 2n scale as S^nM r 1 , characteristic of the shocks. 

A totally different result is found in a system governed by (1) driven by a white-in-time 
random force f(x,t) defined by its correlation function: 



f(k, uo)f(k', u') = 2 (2tt) 2 D k~ y 5(k + k')8(u + w') (2) 

with y = —2 corresponding to thermal equilibrium. Here, in the limit k — > and u — > the 
two-point velocity correlation function is given by 



C(k> a)= n <kMv{kW) ~~~ ~~ <x k-C^) (3) 



with a = z = 3/2 corresponding to E{k)= const. Both exponents z and a were evaluated 
from the theories based on the one-loop renormalized perturbation expansions [||, Q and 
were confirmed by numerical experiments ||. In this case, the small-scale forcing was strong 
enough to prevent formation of the shocks and the fc -2 -energy spectrum. In recent articles 
0, it was shown that computation of the second loops for the y = —2 case does not invalidate 
the results of the one-loop approximation. 

In this work, we are interested in an intermediate case of a system governed by (1) with 
the forcing function added to the right side defined by the relation (2) with y — 1. This 
case is extremely interesting because it corresponds to the "almost" constant energy flux 
H(k) in the wave-number space: H(k) oc log(k/k ), where k Q is the inverse largest allowed 
scale in the system. Since the analytic structure of (1) resembles that of the Navier-Stokes 
equations, the Kolmogorov argument leading to E(k) oc k~s can be applied at least on a 
superficial level. However, in this case the process of the Kolmogorov spectrum generation 
should compete with the natural tendency of the Burgers equation towards formation of 
coherent shocks, thus leading to the most interesting dynamics. 

We investigate the fluctuations generated by the equation (1) with the more general 
dissipation term v$ (— l) p+1 d 2p v/dx 2p and driven by the random force f(x,t). Numerical 
results, shown below, correspond to p — 6, which has been chosen empirically to produce 
sufficiently sharp UV energy spectrum fall-off. The effect of the hyper-dissipation on the 
solutions of the Burgers equation has been studied in a recent paper || and we will not dwell 
upon this issue here. We will just mention that its use is dictated by desire to have as wide 
a universal range as possible and is based on the assumption that universal IR properties 
do not depend on the type of dissipation chosen. To simulate (2), the random force has 
been assigned in the Fourier space as: f(k,t) = Af/y/Si \k\~ y/2 a k , where (jfc is a Gaussian 
random function with | cr fc | =1 and 5t is a time-step. The force cut-off k c was chosen 
well inside the dissipation range of the energy spectrum. In the case p — 1 the dissipation 
scale is, according to the relation given above, ~ jj- where Uq is the rms velocity of the 
most energetic shock in the system. Spatial discretization was based on the Fourier- Galerkin 
pseudospectral method with the nonlinear term computation in the conservative form and 
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dealiasing procedure based on the 1/3-rule. Temporal discretization included two second- 
order schemes: Runge-Kutta for restarting and stiffly-stable Adams-type scheme described 
in || for serial computations. Spectral resolution employed here was 12288 including the 
dealiased modes. Other parameters were chosen to be: u = 9.0 x lCT 40 , 8t = 5.0 x 10" 5 , 
and Af — 1.4142 x 10~ 3 . It was carefully tested that this parameter set does lead to the 
strong coupling in the inertial range 10 < k < 600, such that the viscous term in the energy 
equation derived from (1) is negligibly small compared with the corresponding nonlinear 
term. 

The results of numerical experiments are presented in Figs. 1 — 5. Integration was 
performed for over approximately llr to , where: r to = 7r/V rms ~ 100 is a largest eddy 
turnover time. After approximately 0.5 r to the statistically steady-state has been achieved. 
Fig.l presents 2 successive realizations of the velocity field in this steady-state. One can 
see the typical saw-tooth structures, characteristic of the dynamical system governed by 
Burgers equation. In our case, however, they are superimposed by the random component 
of the velocity field. It was noticed that the system spends most of the time in the state 
where there are few (3-4) large-amplitude shocks and many small-amplitude ones. But 
processes leading to the creation of a single strong shock and its later breakdown into several 
smaller ones constantly take place. The time-evolution of the total energy in the system 
E(t) demonstrates strong (with the amplitude of more than 100% of the average energy) 
fluctuations, characteristic of the instability of the large-scale structures. The time-averaged 
energy spectrum E(k,t) (E(t) = J E(k,t)dk), presented in Fig.2, is well approximated by 
the Kolmoeorov law: E(k) oc k' p with p = § ± 0.02. The error bars were estimated in 
the following way: various values of parameter (3 were used to plot the compensated energy 
spectrum e(k) = k^E(k) and only the values of exponent /3, for which e(k) was inside the 
experimental noise in the entire interval 10 < k < 600, were chosen as represent able. The 
velocity structure functions 5 , 2 n (r) = (v(x + r) — v (x)) 2n are shown in Fig. 3 for n = 2 — 4. 
We can see that all S^r) oc r^ 2n with /5 2n ~ 0.91 indicating that these correlation functions 
are dominated by coherent shocks. We were not able to detect the logarithmic corrections 
to the energy spectrum E(k). However, the fact that the high-order moments, presented 
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on Fig. 3 are characterized by the exponents close but not exactly equal to unity, indicate 
that the logarithmic contributions cannot be ruled out. One remarkable result is related to 



the dissipation rate correlation function presented in Fig.4, G(r)=e(x + r)e(x) oc r~ M , with 
intermittency exponent p, ~ 0.2 ± 0.05 measured inside the universal range 0.01 < r < 0.63. 
Note that the dissipation rate correlation function is defined in the physical space and the 
1/3-rule dealiasing procedure was also used for its computation. The obtained value of the 
dissipation rate exponent p is close to that observed in the experiments on 3-dimensional 
turbulence: p = 0.25 ± 0.05, see flQfl , and its general shape resembles the model shape of 
G(r) proposed for 3-dimensional turbulence in JTIJ. We would like to emphasize that in the 
present work the dissipation rate e(x) = ^(f^) 2 with p = 6 strongly differs from the normal 
viscosity case with p — 1. The fact that the exponent p obtained in this work was close to 
the one observed in real-life turbulence, provides an indication that the correlation function 
of the dissipation rate for the inertial range values of the displacement r is independent 
on the structure of the dissipation range. A similar conclusion was reached from the recent 



numerical experiments of three-dimensional turbulence in ||12|| . As one can see from Fig.4, the 
accuracy of the exponent p is not as good as of the exponent in the expression for the energy 
spectrum. In addition, to assess importance of the result, the role of the hyperviscosity in 
the dissipation rate correlation function is to be further investigated. 

Important information about the dynamics of a non-linear system can be extracted from 
the correlation function defined by (3). In the scale-invariant regime, according to the 
theories based on the one-loop renormalized perturbation expansion [| 

where vik.uj = 0) ~ (^ a )^~3, corresponding to z = 2/3 and a = —7/3 in (3). The 
frequency dependence of the effective viscosity u(k,uj) is neglected in the relation (4). This 
is the direct consequence of the assumption that the dynamics of the inertial range modes 
v(k,u) are dominated by the "distant interactions" with the modes v(q, Q) with \k\ <C \q\ 
and can be decribed by the /c-dependent eddy-viscosity. It is clear that this approximation 
can not be valid when we are interested in behaviour of the most powerfull large-scale 
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structures, because of the strong interaction between them leading to the events of the 
shock instability observed in this work. The energy spectrum derived from (4) is: E(k) = 
J C(k,uo)d uo/(2n) = (ttDq/6) 1 ^ 3 k~ 5 ^ 3 . The energy flux in the wave-number space can be 
expressed in terms of the amplitude of the force correlation function D as folows: n(/c) = 
n(/c ) + 2D log(k/k ). Then, the value of the "Kolmogorov" constant arising from (6) is: 
C K = ((II(fc) -U(k ))/log(k/k )y 2/3 k 5 / 3 E(k) = (tt/6) 1/3 . The numerical value C K « 
0.806 agrees well with the results of numerical simulation, see Fig. 2. As in many other cases, 
understanding of the reasons for a good agreement between the theory, based on the one-loop 
renormalized perturbation expansion, and experimental data, remains a major challenge. 

The computational procedure for C(k,uo) included the following steps. Starting from 
some initial moment t = t well inside the statistically steady-state, the solution v(k,t) 
was stored at the time instants tj = t + T j/M, where T = 100 was chosen to be of the 
order of r to . Then at t — tM v(k,u) was found via discrete Fourier transform. Repeating 
such procedure with time and assuming each realization of v(k,u) to be independent of 
others, which certainly is only an approximation, one can compute C{k,uo) as an average 
over such realizations. Memory requirements allowed us to keep only 200 first wavevectors 
and limit ourselves with M = 3000. Results of computations of C(k,u), presented in Fig.5, 
can be compared with the prediction (4). The relation (4) was derived neglecting the infra- 
red divergences resulting in the transport of the small-scale fluctuations by the large-scale 
coherent structures. This kinematic interaction ( "sweeping effect" ) can be accounted for by 
the Doppler shift uo — > uo + kV in (4), where V is the characteristic velocity of the large 
scale structures. Major interest is whether C(k,uo) is described by (4) or not and whether 
V is zero or not. If the sweeping effect is present in the long-time behavior then there 
are 3 possible scaling regimes of C(k,uo) as uo — > 0+: C(k,u) oc k^ 1 if k (a; 3 j 'Dq) 1 / 2 , 
C(k,uj) oc AT 7 / 3 if (w 3 /£>„) 1/2 < k < D /V 3 , and C(k,u) oc AT 3 if k > D /V 3 . It is 
clear from Fig.5 that the theoretical prediction (4) is surprisingly accurate in both limits 
of large and small frequences u. Only in a narrow intermediate range of the wavenumbers 
uo u(k,uo = 0)k 2 prediction (4) fails. The flattening of C(k,uo) observed in this interval 
indicates that the scaling function F(x) in (3) is a decreasing function of x when x 1. The 
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quantitative agreement between theory and simulations in the limit of large wave-numbers 
k shows that the "sweeping velocity" V is small. This may be a consequence of the fact 
that the large-scale shocks are almost steady. We would also like to note that the accuracy 
of the C(k,u) computation may not be easily increased because of the computer resources 
limits which have been reached by us. The above result leads to an interesting possibility: 
The infra-red divergences, present in the theory, are not summed up into a mere transfer 
of small-scale fluctuations by the large-scale structures but are reflected in the creation of 
a large-scale condensate state, which in this case has a very simple physical meaning as a 
collection of strong shocks moving with a very small velocity V. Derivation of the equation 
of motion describing dynamics of coherent shocks is an important and interesting problem 
and will be the subject of future communication. 

To conclude, the above results, obtained in a simple one-dimensional system are surpris- 
ingly similar to the outcome of experimental investigations of the real-life three-dimensional 
turbulence. In the one-dimensional case, however, the dynamics and geometrical features 
of the "turbulence" building blocks are well understood and a cascade process is readily 
envisioned 3jS du CO 3. gulation of the weak and wide shocks (the shock width and its amplitude 
are related as I ~ ^) into ever stronger and narrower structures until the dissipation takes 
over. Moreover, the total dissipation rate in an interval of the length r is prescribed and 
is equal to e r oc In L ^-. Given these simplifications, one may hope that the full theory of 
Kolmogorov turbulence in the one-dimensional Burgers equation is not out of reach. 

This work was supported by grants from ARPA and AFOSR. Many stimulating discus- 
sions with V. Borue, R. H. Kraichnan, S. Orszag, A. Polyakov and Ya. Sinai are gratefully 
aknowledged. 
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Figure 2: Dotted curve is the energy spectrum E{k) (left axis), the straight line above it has 
the exact slope —5/3. Solid curve is the compensated energy spectrum Ck defined in the 
text (right axis). 
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Figure 3: Velocity differences (v(x + r) — v(x)) 2n for n = 2,3,4 (dotted curve) with the 
linear least-squares fit (solid line). 
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Figure 4: Energy dissipation correlation function e(x + r) e(x) with the linear least-squares 
fit, giving the intermittency exponent \x = 0.2 ± 0.05. 
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Figure 5: Points denote the self-correlation function of the solution C(k,u) for fixed time- 
frequencies uj = lixmjr with m — 1, m — 5, m — 10, m — 15 and r = 100.05. Solid lines 
denote the corresponding asymptotics of the one-loop prediction given by (4). 
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